home *** CD-ROM | disk | FTP | other *** search
/ Aminet 15 / Aminet 15 - Nov 1996.iso / Aminet / dev / basic / ace24dist.lha / ace24.lha / prgs / Astronomy / GalColSim.b next >
Text File  |  1996-09-10  |  14KB  |  570 lines
  1. { Galaxy Collision Simulator 
  2.  
  3.  Galaxy Collision Program
  4.  by M.C. Schroeder and Neil F. Comins.
  5.  Appeared in Astronomy magazine, December 1988
  6.  as an IBM GW-BASIC program. 
  7.  
  8.  Adapted to Macintosh Turbo Pascal 1.1 by David J. Benn
  9.  Date: 4th April - 17th April, 2nd & 3rd May 1991
  10.  
  11.  Adapted for ACE Amiga BASIC by David J. Benn
  12.  Date: 27th,28th November 1992,
  13.        7th January 1993,
  14.        11th October 1993,
  15.        18th December 1993,
  16.        18th September 1994,
  17.        2nd October 1994
  18.     
  19.  This program allows 2 galaxies to interact via
  20.  their gravitational fields. The TARGET galaxy
  21.  is a disk and the INTRUDER galaxy is treated as
  22.  a point mass. }
  23.  
  24. const     SF2=2
  25. const    arr_size=1000
  26. const    aspect=0.625
  27. const    true=-1&
  28. const    false=0&
  29. const    SIM=1
  30. const    PARAMS=2
  31.  
  32. dim     X(arr_size),Y(arr_size),Z(arr_size)
  33. dim    VX(arr_size),VY(arr_size),VZ(arr_size)
  34.  
  35. longint    NRR,NRS,NS,DR
  36.     
  37. single    M1,M2
  38.   
  39. single    X1,Y1,Z1
  40. single    X2,Y2,Z2
  41. single  VX1,VY1,VZ1
  42. single  VX2,VY2,VZ2
  43.     
  44. longint    NTSPR,sim_time
  45.  
  46. SUB title
  47.   title$ = "*** Galaxy Collision Simulator ***"
  48.   color 2
  49.   locate 2,40-len(title$)\2
  50.   PRINT title$   
  51. END SUB
  52.  
  53. SUB await_CR
  54.   color 3
  55.   PRINT
  56.   PRINT "Press <return> to begin simulation..."
  57.   while inkey$<>chr$(13):wend
  58. END SUB
  59.  
  60. SUB instructions
  61. string ans
  62.   PRINT
  63.   color 3
  64.   locate csrlin,25
  65.   PRINT "Do you want instructions? (y/n)"
  66.   repeat
  67.     ans = ucase$(inkey$)
  68.   until ans="Y" or ans="N"
  69.   if ans="Y" then 
  70.     cls
  71.     PRINT
  72.     title
  73.     PRINT
  74.     color 1
  75.     PRINT "This is an adaptation of a BASIC program which appeared in" 
  76.     PRINT "Astronomy magazine, December 1988."
  77.     PRINT
  78.     PRINT "The purpose of the program is to show the large scale effects of"
  79.     PRINT "one galaxy (here refered to as the INTRUDER galaxy) passing nearby"
  80.     PRINT "or through another (TARGET) galaxy."
  81.     PRINT
  82.     PRINT "This is a task which until recent times has been the province of"
  83.     PRINT "of the supercomputer due to the sheer processing power required to"
  84.     PRINT "accurately model the gravitational effects of two galaxies - each"
  85.     PRINT "containing tens or hundreds of billions of stars - upon one ";
  86.     PRINT "another."
  87.     PRINT
  88.     PRINT "However, if the number of stars in the target galaxy is limited and"
  89.     PRINT "the intruder galaxy is considered to be a point gravitational"
  90.     PRINT "source, the problem becomes manageable on a microcomputer."
  91.     PRINT
  92.     PRINT "In this simulation the target galaxy consists of concentric rings" 
  93.     PRINT "of stars."
  94.     await_CR
  95.     cls
  96.     PRINT
  97.     title
  98.     PRINT
  99.     color 1
  100.     PRINT "Before a simulation begins, you will be asked to enter a number of"
  101.     PRINT "parameters. Here is some background information:"
  102.     PRINT
  103.     PRINT "  - 1 parsec = 3.26 light years = approx 30 million million kms"
  104.     PRINT "  - 1 solar mass = the mass of our sun = 1990 thousand million"
  105.     PRINT "    million million million kilograms"
  106.     PRINT "  - Each time step represents 1.2 million years"
  107.     PRINT "  - Each integer X, Y or Z step represents 500 parsecs"
  108.     PRINT "  - Vx, Vy and Vz represent the velocities in the 3 axes"
  109.     PRINT "  - The mass of the target galaxy equals 20 billion solar masses"
  110.     PRINT "  - The mass of the intruder galaxy is entered as a fraction of"
  111.     PRINT "    the target galaxy (eg: 0.25 x 20 billion = 5 billion)."   
  112.     await_CR
  113.     cls
  114.     PRINT
  115.     title
  116.     PRINT
  117.     color 1
  118.     PRINT "The following are some sample initial conditions and their results:"
  119.     PRINT
  120.     PRINT "Result        Intruder Mass   X   Y   Z   Vx   Vy   Vz   Time Steps"
  121.     PRINT
  122.     PRINT "Ring Galaxy       1           7.5  0  35   0    0   -1     65"
  123.     PRINT "Bridge Galaxy     0.25       40   10  10  -1    0    0    120"
  124.     PRINT "Whirlpool Galaxy  0.25      -30   30   0   0   -0.5  0.5  175"
  125.     PRINT
  126.     PRINT "For best results, use about 10 concentric rings and 20 stars per" 
  127.     PRINT "ring. This shows enough detail without slowing the simulation down"
  128.     PRINT "too much."    
  129.     PRINT 
  130.     PRINT "In what follows, you will also be given the opportunity to choose"
  131.     PRINT "one of the above 3 resultant galaxies rather than entering initial"
  132.     PRINT "condition values." 
  133.     PRINT
  134.     PRINT "This program may also be run from the shell/CLI, with the following"
  135.     PRINT "syntax:"
  136.     PRINT
  137.     PRINT "    GalColSim [ring | bridge | whirlpool]"
  138.     await_CR
  139.   end if
  140. END SUB
  141.  
  142. SUB entry_method
  143. string ans 
  144.   cls
  145.   color 1
  146.   PRINT:PRINT
  147.   PRINT "Enter [P]arameters or [N]ame of simulation?"
  148.   repeat
  149.     ans=ucase$(inkey$)
  150.   until ans="P" or ans="N"
  151.   case
  152.     ans="P" : entry_method=PARAMS
  153.     ans="N" : entry_method=SIM
  154.   end case
  155. END SUB
  156.  
  157. SUB set_parameters(opt)
  158. shared      NRR,NRS,NS,DR
  159. shared     M2
  160. shared     X2,Y2,Z2
  161. shared     VX2,VY2,VZ2
  162. shared     NTSPR
  163.   '..initialise parameters
  164.   if opt=1 then
  165.     '..ring
  166.     M2=1
  167.     X2=7.5 : Y2=0  : Z2=35
  168.     VX2=0  : VY2=0 : VZ2=-1    
  169.     NTSPR=65
  170.   else
  171.     if opt=2 then
  172.       '..bridge
  173.       M2=0.25
  174.       X2=40  : Y2=10 : Z2=10
  175.       VX2=-1 : VY2=0 : VZ2=0
  176.       NTSPR=120   
  177.     else
  178.       '..whirlpool
  179.       M2=0.25
  180.       X2=-30 : Y2=30    : Z2=0
  181.       VX2=0  : VY2=-0.5 : VZ2=0.5
  182.       NTSPR=175
  183.     end if
  184.   end if 
  185.  
  186.   '..number of rings and stars per ring
  187.   NRR=10
  188.   NRS=20
  189.   NS=NRR*NRS
  190.   if NRR<2 then NRR=2
  191.   DR=20 \ (NRR-1)
  192. END SUB
  193.  
  194. SUB args_present
  195. string end_state
  196.   if argcount<>1 then
  197.     args_present=false
  198.     exit sub
  199.   else 
  200.     args_present=true
  201.     end_state=ucase$(arg$(1))
  202.     case
  203.       end_state="RING"      : opt=1
  204.       end_state="BRIDGE"    : opt=2
  205.       end_state="WHIRLPOOL" : opt=3
  206.     end case
  207.     set_parameters(opt)       
  208.   end if
  209. END SUB
  210.  
  211. SUB get_simulation
  212. shortint opt
  213.  
  214.   '..ask for simulation end-state
  215.   PRINT
  216.   PRINT "Which end-state?"
  217.   PRINT
  218.   PRINT "1. Ring Galaxy"
  219.   PRINT "2. Bridge Galaxy"
  220.   PRINT "3. Whirlpool Galaxy"
  221.   repeat 
  222.     opt=val(inkey$)
  223.   until opt>=1 and opt<=3
  224.   
  225.   set_parameters(opt)
  226. END SUB
  227.  
  228. SUB get_parameters
  229. {The stars are distributed in the
  230.  TARGET galaxy in rings. The total number
  231.  of stars is the product of the # of rings and
  232.  the # of stars per ring.}
  233. shared     NRR,NRS,NS,DR
  234. shared    M2
  235. shared    X2,Y2,Z2
  236. shared    VX2,VY2,VZ2
  237. shared    NTSPR
  238. string     ANS  
  239.  repeat
  240.   cls
  241.   PRINT
  242.   PRINT "Enter the initial conditions..."
  243.   PRINT
  244.   PRINT "Input the number of rings of stars in the TARGET galaxy (try 8-20):"
  245.   input NRR
  246.   PRINT
  247.   PRINT "Input the number of stars per ring in the TARGET galaxy (try 15-30):"
  248.   input NRS
  249.   PRINT
  250.   NS=NRR*NRS
  251.   if NRR<2 then NRR=2
  252.   DR=20 \ (NRR-1)
  253.   PRINT "Input the mass fraction of the INTRUDER galaxy"
  254.   PRINT "in units of the TARGET galaxy mass (try 0.25-1):"
  255.   input M2
  256.   PRINT
  257.   PRINT "Input the initial X Y Z co-ordinates of the INTRUDER galaxy"
  258.   PRINT "(eg: 40 10 10):"
  259.   PRINT "The TARGET galaxy is located at 0 0 0."
  260.   PRINT "X: "; 
  261.   input X2
  262.   PRINT "Y: "; 
  263.   input Y2
  264.   PRINT "Z: "; 
  265.   input Z2
  266.   PRINT
  267.   PRINT "Input the initial X Y Z  velocities"
  268.   PRINT "of the INTRUDER galaxy (eg: -1 0 0):"
  269.   PRINT "The TARGET galaxy is initially at rest."
  270.   PRINT "X: "; 
  271.   input VX2
  272.   PRINT "Y: "; 
  273.   input VY2
  274.   PRINT "Z: "; 
  275.   input VZ2
  276.   {not too small (0 is ok) - system may crash at around +- .34 or so...}
  277.   if VX2<0.5 then VX2=VX2*2;
  278.   if VY2<0.5 then VY2=VY2*2;
  279.   if VZ2<0.5 then VZ2=VZ2*2;
  280.   PRINT
  281.   {the position of the TARGET galaxy
  282.    is assigned by the computer}
  283.   PRINT "Input number of time steps for this run (try 50-200)"
  284.   input NTSPR
  285.   PRINT
  286.   PRINT "Are these inputs correct? (y/n) "
  287.   repeat
  288.     ANS = ucase$(inkey$)
  289.   until ANS="Y" or ANS="N"
  290.  until ANS="Y"
  291. END SUB
  292.     
  293. SUB update_display
  294. {update screen display}
  295. shared     X,Y,VX,VY,Z,VZ
  296. shared    NRR,NRS,NS,DR    
  297. shared    M1,M2  
  298. shared    X1,Y1,Z1
  299. shared    X2,Y2,Z2
  300. shared  VX1,VY1,VZ1
  301. shared  VX2,VY2,VZ2   
  302. shared    NTSPR,sim_time
  303. single     XC,YC,ZC,XX1,YY1,ZZ1,XX2,YY2,ZZ2
  304. longint    I,J
  305. single    XP,YP,ZP
  306. string    ch
  307.     {calculate centre of mass of system
  308.      for use as centre of output display}
  309.     XC=(M1*X1+M2*X2)/(M1+M2)
  310.     YC=(M1*Y1+M2*Y2)/(M1+M2)
  311.     ZC=(M1*Z1+M2*Z2)/(M1+M2)
  312.     {set up output display - clear text & graphics}
  313.     cls
  314.     {print out display headings}  
  315.     color 3
  316.     locate 2  
  317.     PRINT "        Polar             Edge-on"
  318.     color 3
  319.     locate 22
  320.     PRINT "        Step:";sim_time
  321.     color 5
  322.     locate 23
  323.     PRINT "   Real Time:";int((sim_time+1)*1.2);"million years"
  324.     color 4
  325.     locate 24
  326.     PRINT "   Press any key to stop"
  327.     '..view frames
  328.     LINE (0,0)-(159,166),2,b
  329.     LINE (161,0)-(319,166),2,b
  330.     {calculate position of galactic centres and display on screen}
  331.     XX1=(X1-XC)
  332.     YY1=(Y1-YC)
  333.     ZZ1=(Z1-ZC)
  334.     XX2=(X2-XC)
  335.     YY2=(Y2-YC)
  336.     ZZ2=(Z2-ZC)
  337.     ZP=Z(1)-ZC + Z(NS)-ZC
  338.     CIRCLE(int(XX1+75),int(YY1+75)),2,1,,,aspect    {target x-y}
  339.     CIRCLE(int(XX1+225),int(ZP+75)),2,1,,,aspect    {target x-z}
  340.     CIRCLE(int(XX2+75),int(YY2+75)),1,5,,,aspect         {intruder x-y}
  341.     CIRCLE(int(XX2+225),int(2*ZZ2+75)),1,5,,,aspect       {intruder x-z}
  342.     {place stars on screen}
  343.     FOR I=1 TO NS
  344.       XP=int(X(I)-XC)
  345.       YP=int(Y(I)-YC)
  346.       ZP=2*(Z(I)-ZC)
  347.       PSET(int(XP+75),int(YP+75)),1
  348.       PSET(int(XP+225),int(ZP+75)),1
  349.     NEXT
  350. END SUB
  351.  
  352. SUB define_galaxies
  353. shared     X,Y,VX,VY,Z,VZ
  354. shared    NRR,NRS,NS,DR    
  355. shared    M1,M2  
  356. shared    X1,Y1,Z1
  357. shared    X2,Y2,Z2
  358. shared  VX1,VY1,VZ1
  359. shared  VX2,VY2,VZ2   
  360. shared    NTSPR,sim_time
  361. longint    IR,IT,I
  362. single    R,RV,T,T1,TH,V
  363.  PRINT
  364.  PRINT "*** Defining structure & position" 
  365.  PRINT "    of galaxies. Please wait."
  366.  {initialise TARGET galaxy mass,
  367.   position and velocity}
  368.  M1=5
  369.  X1=150
  370.  Y1=100
  371.  Z1=0
  372.  VX1=0
  373.  VY1=0
  374.  VZ1=0
  375.  {scale the INTRUDER galaxy mass,
  376.   position and velocities.}
  377.  sim_time=0
  378.  M2=M2*M1
  379.  X2=X2+X1
  380.  Y2=Y2+Y1
  381.  Z2=Z2+Z1
  382.  {setup initial load
  383.   place stars in concentric rings
  384.   from r=10 to r=30 at intervals of DR}
  385.  FOR IR=1 TO NRR
  386.     R=10+(IR-1)*DR
  387.     V=SQR(M1/R)
  388.     TH=(0.5*V/R)*(180/3.14159)
  389.     IF R=10 THEN V=0.9*V
  390.     FOR IT=1 TO NRS
  391.         T=(IT-1)*360/NRS
  392.         T1=3.14159*(T-TH)/180
  393.         I=I+1      
  394.         {initialise star positions}
  395.         X(I)=R*COS(T/57.2958)+150
  396.         Y(I)=R*SIN(T/57.2958)+100
  397.         VZ(I)=0
  398.         Z(I)=0
  399.         {initialise star velocities to place them in circular
  400.          orbits about the TARGET galaxy.}
  401.         VX(I)=-V*SIN(T1)
  402.         VY(I)=V*COS(T1)
  403.     NEXT
  404.  NEXT
  405. END SUB
  406.  
  407. SUB particle_pusher
  408. shared     X,Y,VX,VY,Z,VZ
  409. shared    NRR,NRS,NS,DR    
  410. shared    M1,M2  
  411. shared    X1,Y1,Z1
  412. shared    X2,Y2,Z2
  413. shared  VX1,VY1,VZ1
  414. shared  VX2,VY2,VZ2   
  415. shared    NTSPR,sim_time
  416. longint    I,J,K
  417. single    R1,R2
  418. single    AX,AY,AZ
  419. single    RR 
  420. single    Ex,Ey,Ez
  421.     FOR K=1 TO NTSPR 
  422.         FOR J=1 TO 1 
  423.             FOR I=1 TO NS
  424.              {determine the force on a star
  425.               due to the galactic centres
  426.               SF2 below, called the softening factor,
  427.               is used to prevent overflows
  428.               during force calculations
  429.               SF2 is assigned a value above.
  430.               You may experiment with its value
  431.               but it should be kept as small as possible
  432.               to better approximate a true 1/r**2 force.}
  433.               Ex=(X(I)-X1)^2
  434.               Ey=(Y(i)-Y1)^2
  435.               Ez=(Z(i)-Z1)^2
  436.               R1=M1/(Ex+Ey+Ez+SF2)^1.5
  437.               Ex=(X(I)-X2)^2
  438.               Ey=(Y(i)-Y2)^2
  439.               Ez=(Z(i)-Z2)^2
  440.               R2=M2/(Ex+Ey+Ez+SF2)^1.5
  441.               {calculate stars x,y,z accelerations}
  442.               AX=R1*(X1-X(I))+R2*(X2-X(I))
  443.               AY=R1*(Y1-Y(I))+R2*(Y2-Y(I))
  444.               AZ=R1*(Z1-Z(I))+R2*(Z2-Z(I))
  445.               {update star positions and velocities
  446.                using a time centred leap-frog algorithm}
  447.               VX(I)=VX(I)+AX
  448.               VY(I)=VY(I)+AY
  449.               VZ(I)=VZ(I)+AZ
  450.               X(I)=X(I)+VX(I)
  451.               Y(I)=Y(I)+VY(I)
  452.               Z(I)=Z(I)+VZ(I)
  453.               if inkey$<>"" then exit sub
  454.             NEXT
  455.           {update positions
  456.            and velocities of the galactic centres}
  457.            Ex=(X1-X2)^2
  458.            Ey=(Y1-Y2)^2
  459.            Ez=(Z1-Z2)^2
  460.            RR=(Ex+Ey+Ez+SF2)^1.5
  461.            AX=(X2-X1)/RR
  462.            AY=(Y2-Y1)/RR
  463.            AZ=(Z2-Z1)/RR
  464.            VX1=VX1+M2*AX
  465.            VY1=VY1+M2*AY
  466.            VZ1=VZ1+M2*AZ
  467.            VX2=VX2-M1*AX
  468.            VY2=VY2-M1*AY
  469.            VZ2=VZ2-M1*AZ
  470.            X1=X1+VX1
  471.            Y1=Y1+VY1
  472.            Z1=Z1+VZ1
  473.            X2=X2+VX2
  474.            Y2=Y2+VY2
  475.            Z2=Z2+VZ2
  476.            ++sim_time
  477.         NEXT
  478.         update_display
  479.      NEXT   
  480. END SUB
  481.  
  482. SUB continue
  483. shared     NTSPR
  484. string     ans
  485.     color 2
  486.     locate 24
  487.     PRINT "   Continue simulation? (y/n)"
  488.     repeat
  489.       ans = ucase$(inkey$)
  490.     until ans="Y" or ans="N"
  491.     if ans="Y" then 
  492.       continue=true
  493.       locate 24
  494.       PRINT "   Enter number of time steps: ";
  495.       input NTSPR
  496.     else 
  497.       continue=false
  498.     end if
  499. END SUB
  500.  
  501. SUB finished 
  502. string ans
  503.   cls
  504.   color 5
  505.   locate 12,8
  506.   PRINT "Another Simulation? (y/n)"
  507.   repeat
  508.     ans=ucase$(inkey$)
  509.   until ans="Y" or ans="N"
  510.   if ans="N" then finished=true else finished=false
  511. END SUB
  512.  
  513. SUB open_hires
  514.   SCREEN 2,640,225,2,2
  515.   WINDOW 2,,(0,0)-(640,225),32,2
  516.   palette 0,0,0,0    '..black
  517.   palette 1,1,1,1    '..white
  518.   palette 2,1,1,0.13    '..yellow
  519.   palette 3,1,0,0    '..red
  520.   FONT "topaz",8
  521.   color 1
  522. END SUB
  523.  
  524. SUB close_hires
  525.   WINDOW CLOSE 2
  526.   SCREEN CLOSE 2
  527. END SUB
  528.  
  529. { ** main ** }  
  530. SCREEN 1,320,225,3,1
  531. WINDOW 1,,(0,0)-(320,225),32,1
  532. palette 0,0,0,0        '..black
  533. palette 1,1,1,1        '..white
  534. palette 2,1,0,0        '..red
  535. palette 3,0,1,0        '..green
  536. palette 4,0,0,1        '..blue
  537. palette 5,1,1,0.13    '..yellow
  538.  
  539. FONT "topaz",8
  540.  
  541. color 1
  542.  
  543.     argument = args_present
  544.  
  545.     repeat  
  546.       if not argument then 
  547.         open_hires
  548.         title
  549.         instructions 
  550.         if entry_method = SIM then
  551.           get_simulation
  552.         else
  553.           get_parameters
  554.         end if
  555.         define_galaxies 
  556.         close_hires
  557.         update_display
  558.       else
  559.         define_galaxies
  560.     update_display
  561.         argument=false
  562.       end if
  563.       repeat   
  564.         particle_pusher
  565.       until not continue
  566.     until finished
  567.  
  568. WINDOW CLOSE 1
  569. SCREEN CLOSE 1
  570.